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Abstract 

Using the archives of the American Association of Variable Stars Observers 
and our own data, we analyse the long-term variability of several well-studied 
Luminous Blue Variables (LBVs) aiming on a general picture of stochastic 
variability of these objects. The power density spectra of all the selected 
objects may be generally described by a single power law contaminated by 
observational noise at higher frequencies. The slopes of the power-law com- 
ponent are close to p = 2 (where PDS oc /^^, and / is frequency) for strongly 
variable flaring objects like AG Car and significantly smaller (p ~ 1.3) for 
P Cyg where brightness variation amplitude is < 1™ and dominated by slow 
low-amplitude variability. The slope holds for about two orders of magni- 
tude in the frequency domain, though peaks and curvatures are present at 
/ ~ 10~^-i-10~'^d~^. We show that pseudo-photosphere approach to variabil- 
ity may explain the power-law shape of the variability spectrum at higher 
frequencies. However, the observed spectra are actually rather "red" than 
"brown" : flux variations are correlated up to tens of years that is much longer 
than the characteristic refreshment time scales of the pseudo-photosphere. 
We propose that several stochastic noise components produce the power spec- 
tra of LBVs. 

Keywords: stars: activity, stars: winds, outflows, stars: individual (P 
Cygni, eta Carinae, AG Carinae, Romano star, M33 V* V268 (=var C) ) 



1. Introduction 



Flares from luminous blue variables (LBV) are studied for about four hun- 
dred years, since the 1600 flare of P Cyg (de Groot et al. , 2001). However, 



these objects are rare and thus relatively unstudied. As variable objects, 
LBV are characterised by a complex hierarchy of variability timescales and 
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amplitudes (Humphreys and Davidson 1994). Usually, three empirical vari- 
ability scales are distinguished: 



1. microvariations, with amplitudes ~ 0™1 and relatively small variability 
timescales, from days to months 

2. S Dor cycles, or eruptions (following Humphreys and Davidson (1994)) 
with amplitudes ~ 1 -i- 2^, years to tens of years in length 

3. giant eruptions with similar and larger (up to ~ lOOyr) timescales and 
larger amplitudes (> 2^) 

Physically, giant eruptions may differ from the variations of the second 
type by significant changes in the bolometric luminosity (see discussion in 
section 5.2). Some authors like van Genderen et al. ( 1997a[ ) propose division 
of the second type variability into normal and very long timescale S Dor 
cycles. Longer-timescale variations indeed seem to represent a separate mode 
of LBV variability, that is most evident in the XXth-century light curve of 
P Cyg (see next section) that is free from ordinary fiares but dominated by 
microvariations superimposed over slow variations > lOyr in length. 

Several well-studied LBVs exhibit periodic variability components. In 
particular, van Genderen et al. (1997b) mention a ~ lyr period for AG Car 
and a ~ 7yr period for S Dor. Probably, these periods do not correspond to 
any coherent process but rather refiect some characteristic timescales close 
to the duration timescale of a single fiare. 

7] Carinae is an exceptional case in this sense. The well-known 5-year 
period is stable on the timescales of tens of years that supports the idea 



that this object is indeed a (rather broad) binary (Damineli et al. , 2008). 
Note however that the orbital separation of a massive binary with a 5-year 
period should be a ~ 2 x 10^^{M/lOOMQY/^iT /byrY/^cm. Hydrostatic radii 
for LBV and related stars are ~ lO^^cm that makes the expected effects of 
binary interaction on the central machine (-s) of rj Car very small. In this 
study, we find the broad-band PDS of rj Car quite similar to that of other 
fiaring LBVs such as AG Car. 

All the observed variability time scales of LBV stars lie between their 
dynamic and thermal times. Dynamic time scale is determined by the mean 
density of the stellar hydrostatic core. Outflowing matter can not take part 
in pulsations because the wind rapidly becomes supersonic. Straightforward 
estimate for the dynamic timescale as the free-fall time yields: 
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Due to period doubling (see for example Buchler (1993)) and instability 
of pulsational modes, lower-frequency modes are excited during pulsations. 
Probably, this effect is observed in "ordinary" hot supergiant stars, where 



microvariations occur at unexpectedly long time scales of several days (van 



Genderen 1985). Flaring LBV variability is still much longer. 

On the other hand, Kelvin-Helmholtz time scale for a typical hypergiant 

is: 

Unfortunately, these time scales are unreachable for modern observations, 
but may be, in principle, studied indirectly by observations of the ejecta, 
possibly light echoes from strong flares and using the statistical properties of 
LBV ensembles. 

All the considered variability scales are somewhere in between and may 
correspond either to secular effects in dynamic-timescale evolution ("slow 
pulsations") or to effects in the expanding atmosphere that has, evidently, 
longer characteristic timescales than the hydrostatic core. Stellar wind, as- 
cribed some velocity Vw, mass loss rate M and opacity /t, is characterised by 
the radius of the (pseudo-) photosphere Rph — i^M /AttVw and corresponding 
time scale that has the physical meaning of the replenishment time for the 
matter in the wind: 

W - R.^h^ - 2;!^10-5Meyr- ( 100^) ^ ^ 

Here, kt — 0.4cm^g~^ is the Thomson cross-section for the scattering by 
free electrons in a hydrogen-rich gas. Real opacity is higher and contains some 
contribution of true absorption processes that may increase the characteristic 



time scales by about an order of magnitude (see, for example, Iglesias and 



Rogers (1996), bearing in mind that, for a typical LBV, IgT ~ 4 5 and 
Igi? ~ —5 -. — 7 for the inner wind). Despite the apparently small value 
of tj-gjj it is an important estimate for the characteristic variability time of 
an optically thick hot stellar wind. Longer time scales arise due to higher 
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opacities Mass loss modulation during strong outbursts leads to even longer 
time scales, up to years. 

Bolometric luminosity is not conserved during strong flares (see discussion 
in section 5.2). Besides this, modeling LBV photospheres provides evidence 



(de Koter et al. , 1996) for more complicated nature of the variability of these 



objects, incorporating pulsations, mass loss variations and some additional 
mechanisms (such as unstable pulsational modes) leading to tremendous en- 
ergy release during giant eruptions. 

The overall power density spectrum (PDS) of LBV variability has not 
been considered as a whole. It is tempting to compare LBV stars to active 
galactic nuclei where different time scales of a single flicker noise were for 
a long time interpreted as physically distinct variability processes (see for 



example Terebizh et al. (1989)). Below we show that indeed for several LBV 



stars, broad-band power spectra have steep power-law shapes. 

In this work we aim for the broad-band PDS of LBV variability. In 
the next section, we make a compilation of observational data on several 
reasonably well-studied LBV stars and present their PDS in the following 
section |3| A possible explanation for the observed broad-band PDS shapes 
is proposed in section |4j Results are discussed in section [5| 



2. Observational Data 

We selected three well-studied and representative Galactic LBV stars. For 
these objects (see table [T]), relatively long observational series exist, spanning 
periods of time in excess of 50 years. The data were taken from the archives 
of the American Association for Variable Star Observers. Earlier data are 
primarily visual magnitude estimates, therefore we use only visual magni- 
tudes. We also checked that the data series do not have gaps longer than 
several months and binned the observational data points by five to diminish 
the effect of single erroneous estimates and the round-off effect connected 
to the low precision (~ CPl) of visual magnitude measurements. Note that 
using a median filter would not diminish the round-off effect. 

We also consider V-band light curves of comparable length obtained for 
two Hubble- Sandage variables in M33, Romano's star (the light curve itself 



(Humphreys et al. , 1988) 



is described in Sholukhova et al. (2002)) and Hubble-Sandage Variable C 
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Table 1: Objects selected for analysis. For Galactic objects, the numbers of data points 
are given after binning. 



Object ID 



Variability Limits, mag 



Time span 



Number of Data Points 



P Cyg 
r] Car 
AG Car 

Romano's star 
Var C 



Galactic objects 

4.5^ 5.6 1917-08-13.. 2010-05-07 

4.5^ 8.0 1943-07-24.. 2010-07-19 

5.6H- 8.9 1939-12-05.. 2010-05-03 

M33 

16.1^ 18.6 1949-12-20. .2009-11-08 

15.3^ 17.9 1961-09-13.. 2005-11-08 



3396 
4637 
2186 

802 
635 



2001). 



2.1. Notes on Individual Objects 
P Cyg 

P Cyg is known for its outbursts in 1600 and later (de Groot et al. 
However, in the considered nearly 100-year time span it hardly shows any 
signatures of activity, except for rare short-lasting low-amplitude excursions 
(Markova et al. , 2001). Variability amplitude is primarily contributed by 
low-frequency variability. Several periodicities, from 17 to 100 days, were 
reported (iKolkal 120011). 



rj Car 

During the considered time span, its behaviour is as well relatively quiet, 
especially if compared with its tremendous outbursts in the XIX century. In 
our sample, r] Car is the only proven binary. Its orbital period does not show 
up strongly in the optical variability, probably because of the large optical 
depth of the outflows in the optical. The variability pattern is dominated by 
the gradual increase of the optical luminosity. 



AG Car 



Variability of this prototypical LBV star was considered in [van Genderen 



et al. (1997c) together with that of r] Car. The object exhibits several strong 
flares. Periodicity of ~ 370^^ (somewhat smaller than the characteristic flare 
duration time scale) is known to be present. 
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Romano's Star, or V532 (M33) 

One strong flare with an amplitude of about 2" and a couple of smaller 
maxima represent a typical picture of S Dor variability. The V-band data and 



data reduction process are described in Sholukhova et al. (2002). The object 



is hotter than all the other stars of the sample and reaches the spectral class 



of WN8 in its low/hot state (Polcaro et al. , 2011; Maryeva and Abolmasov 



2010). Its spectrum in the quiet state (Bla^e) is similar to that of P Cyg. 



Var C (M33) 



Behaviour of this variable denoted as Variable C by Hubble and Sandage 



(1953) was studied by Humphreys et al. (1988), who conclude that the object 



behaves as a typical S Dor variable with a higher than average mass loss rate. 
The variable spends about a half of its time in a high and cool phase of the 
S Dor cycle. Amplitudes of the flares reach 2™5 in the visual band. The data 
were provided by Alia Zharova (private communication). 



3. Analysis and Results 

All the time series are non- uniform but free from strong aliases, therefore 



we use extirpolation method (Press and Rybicki, 1989), where the initially 



non-uniform time series is interpolated upon a regular grid, to compute the 
Fourier spectra. Fourier transform for uniform series was performed using 
the fftw library (Frigo and Johnson, 2005). Software written in C and IDL 



was used for extirpolation and binning of the power density spectra (PDS). 
PDS were binned by 10. We use relative normalization defined as follows: 



1 



AT 



(4) 



where {F) is the mean flux, Fj is the j-th component of digital Fourier 
transform of the flux F (already interpolated over a regular grid), AT = iVAi 
is the total effective time span, N is the number of data points, j = 0..N — 1, 
At is the spacing of the regular grid. This normalisation has an evident 
physical meaning of the relative variability power in a unit frequency range. 
Its expectation does not depend on the time span considered and on data 
binning. Variability amplitude may be estimated by integrating the PDS 
over frequency domain and taking a square root. 
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Table 2: PDS Fitting Results. 
Galactic 

P Cyg rj Car AG Car 



M33 

V532 Var C 



TV, 10-^-1 1.9 ±0.2 

a 1.17±0.02 
X^/DOF 549/168 



Power Law: 

5.3 ±0.4 1.6 ±0.3 
1.31 ±0.01 1.56 ±0.04 
1393/230 545/107 



0.59 ± 0.17 0.42 ±0.09 
1.63 ± 0.05 1.86 ±0.05 
166/38 213/30 



Lorentzian: 



TV, d-i 
7, 10-M-i 
X^/DOF 



38 ±4 
2.16 ±0.07 
1602/168 



770 ± 30 
1.69 ±0.03 

2574/230 



1900 ± 1200 
1.4 ± 1.3 
2167/107 



1.23 ± 1.07 
1.23 ± 1.07 
169/38 



1.51 ±0.18 
1.51 ±0.18 
227/30 



Power Law with White Noise Component: 

TV, 10-M-i 0.47 ± 0.07 0.67 ±0.06 0.090 ±0.007 0.008 ± 0.004 0.08 ± 0.03 

a 1.43 ± 0.03 1.75 ±0.02 2.12 ±0.01 2.34 ±0.01 2.14 ±0.07 

TVvF, 10"*d-i 45.2 ±4.9 139 ± 6 302 ± 15 550 ± 20 500 ± 100 

X^/DOF 536/167 1344/229 455/106 123/37 199/29 



All the obtained PDS (see figure [T]) have approximately power-law shapes 
with putative flattenings at lower (< 10~^d~^) and higher (> O.Old"^) fre- 
quencies. The lower-frequency turnovers are smoothed out by the binning in 
the frequency domain, but are better seen at finer binning. The data were fit- 
ted with power-law {PDS = A^- /""), Lorentzian {PDS = N / {l + {fhf), 
with zero resonance frequency) and power-law -|- white noise {PDS = N ■ 
/~" + Nq) models. Our choice of the third model is justified by existence of 
the observational errors that create an uncorrelated (white) noise component. 
Lorentzian function is relevant because it is expected to represent the PDS of 
a Poissonian sequence of exponentially decaying flares (and, approximately, 
rapidly-decaying flares of other shapes; see section 4.4). It indeed proves to 
be a reasonable fit for the objects that demonstrate strong flaring activity. 
Fitting results are given is table [2j 

P Cyg here is unique in its spectral shape that has the hardest power- 
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Figure 1: Power density spectra of the five variables fitted by the power-law + white noise 
model. 
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law slope spanning a frequency range of two orders of magnitude. It is also 
unique in being the least active among all the objects. Other stars that 
exhibit S Dor cycles and outbursts have steeper spectral slopes close to 2 
with possible indications for a low-frequency flattening. Average LBV power 
spectrum clearly has a power-law shape between tens of days to years or tens 
of years and possibly further toward hundreds of years. Gaining statistics and 
increasing the time spans to hundreds of years may allow to clearly resolve 
the turnover at about / ~ 10~^ d~^ that definitely should take place because 
the light curve should remain reasonably uniform in time. 

All the spectra exhibit curvatures and diffuse peaks that make all the fits 
rather poor {x^ / DOF ~ 2-f-5). The overall broad-band shape is, however, in 
good agreement with a single power law (with a high-frequency contribution 
from observational uncertainties) in all the cases. 



4. Variability Pattern of a Pseudo-Photosphere 

The idea of pseudo-photospheric nature of LBV cycles may be traced 



down to the note of Lamers (1987) who proposed that the primary source of 
the variability in these objects is the variability of the mass-loss rate. Below, 
we will assume that this variability is fast and uncorrelated (white noise in 
mass loss rate) and show that such variations indeed lead to a Brownian- 
like noise in the PDS, but with a turnover at a higher frequency. Lower- 
frequency parts of LBV PDS require some additional mechanism of stochastic 
variability. 

^.1. Basic Assumptions 

Excluding brightest outbursts, variability of luminous blue variables is 



generally assumed to be purely spectral. According to Lamers (1987), the 
basic mechanism responsible for variable luminosities of LBV stars in the 
optical is through variations in the mass loss rate resulting in variable pseudo- 
photosphere radius and effective temperature. 

In the relatively hot, relatively rarefied atmospheres of massive stars opac- 
ities are dominated by Thomson scattering, therefore the (radial) optical 
depth through the wind may be calculated as: 

.(j;.,)^^r^^"7/t»* (6) 

47r Jji r''v{r) 
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Here, v = v{r) is the velocity of the wind. Different true absorption 
processes also contribute to the opacity of the wind, but their effect depends 
on the temperature structure of the wind and on the effects of clumping. We 
neglect all the effects of true absorption though they are able to increase the 
opacity of the wind by a factor of several. This will mimic a larger mass 
accretion rate in the model. In its simplest form, the radius and effective 
temperature as a function of time may be estimated by equating the optical 
depth to some fixed value (below we use r = 2/3). Equation ^ may be then 
solved for r. 



4-2. Numerical Simulations 

Six model light curves (see table [s]) were calculated using a pseudo- 
photosphere model with a blackbody photosphere defined by the condition of 
T = 2/3. Bolometric luminosity was everywhere set to L^oi = 3 x lO^^ergs"^ 
that is close to its value for P Cyg as estimated by Pauldrach and Puis (1990). 
All the curves are 50 years in length and contain between 195 and 9774 data 
points. Model F was aimed on reproducing the basic variability properties 
of P Cyg. Others differ in the unperturbed mass loss rate Mq, mass-loss rate 
dispersion and wind velocity. 

We use the semi-empirical /3-law v{r) = {voo — fo) (1 — r/R)^ + vq with 
(3 = 1. Wind acceleration spatial scale R is assumed identical to the hydro- 
static inner radius R, and we set vq = 2.5km s~^ everywhere. Wind velocity 
at infinity was fixed, while the mass loss rate experiences log-normal varia- 
tions that we assume uncorrelated on the time scales of interest. Of course, 
real stellar winds have variable velocities, but simultaneous variations of wind 
velocity and mass loss rate do not change the picture qualitatively. Variable 
wind velocity would however lead to formally divergent solutions for density 
distribution. In real winds, shock waves are expected to be formed. Hence 
the effects of variable wind velocity should be considered together with de- 
viations from spherical symmetry and quasi-stationarity. Log-normality of 
the mass-loss rate distribution is a useful assumption because is allows to 
easily account for dramatic changes in the value of M without producing 
unphysical negative values. If the mass-loss rate varies on dynamical time 
scales, the wind will work as an integrator, making the observed power den- 
sity spectrum softer than the spectrum of M variations. For most of the 
models, we assume the mass loss rate logarithm dispersion a (In rh) = 2 that 
corresponds to about a factor of 7 change in the mass loss rate itself and may 
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Table 3: Light Curve Simulations 



Model ID 


A 


B 


C 


D 


E 


F 


Mo, Moyr-i 


10-4 


10-4 


10-5 


10-4 


10-4 


10-5 


a (Inm) 


2 


2 


2 


3.5 


1 


1 


■y / 100km 


1 


0.5 


1 


1 


2 


2 


i?o, IQi^cm 


28 


14 


3.7 


28 


55 


3.7 


No of points 


3910 


1955 


9774 


195 


1955 


5585 


My range, mag 


-6.6..-10.0 


-6.7..-10.0 


-5.7..-8.2 


-6.1. .-10.0 


-6.6. .-10.0 


-5.7..-7.9 


(My), mag 


-8.2 


-9.0 


-6.0 


-9.0 


-9.0 


-6.0 


cr(My), mag 


0.5 


0.5 


0.24 


0.67 


0.5 


0.22 


N{tdyn < tdiff), % 


0.2 








4 


0.2 






account for the observed ~ 1 2" variations of S Dor cycles. Models D and 
E illustrate the effect of variable amplitude of mass loss variations. 

Inner (hydrostatic) radius was everywhere set to 2 x lO^^cm, the scale 
height of the /3 law for velocity is ascribed the same value. 

In figure [2| we show the simulated light curves in the Johnson V band. 
Corresponding power density spectra are shown in figures [3] and |4] together 
with modified power law and Lorentzian approximations. Modified power 
law is defined as follows: 



'MPL 



Nrexp{-Tf) + No (6) 



This spectral law differs from the power law + white noise model used 
above by an exponential decay factor with characteristic time T, close to 
the wind replenishment time. If the mass loss is uncorrelated at higher 
frequencies, small fast variations should have a flat spectrum as well (see 
section 4.3). 



Basic parameters of the light curves and fitting results are given in tables 
|3]and|4| respectively. Nit^yn < tdif) in the table is the fraction of time when 
dynamical timescale defined as tdyn = Ro/v is smaller than the radial diffu- 
sion time scale defined as t^iff = tRq/c, where r is the total optical depth 
throughout the wind. If this inequation is satisfied, the assumption of in- 
stantaneous energy transfer vital for the pseudo-photosphere approximation 
is violated (see section |5.2[ ). Power density spectra of the simulated light 
curves are shown in figures [3] and |4] 
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Figure 2: Simulated light curves. 



Table 4: PDS of the Simulated Light Curves 



Model ID 



B 



C 



D 



E 



Lorentzian: 



N, d 

7, lO-^d" 
XVDOF 



30 ± 1 
1.44 ±0.05 
2527/115 



0.468 ± 0.007 

3.92 ±0.19 
1304/298 



0.91 ±0.09 

14 ±3 
1299/291 



120 ± 10 
1.2 ±0.3 
55/24 



42 ±2 
0.9 ±0.1 
947/27 



1.2 ±0.1 

6.48 ± 0.05 
1700/333 



Modified Power Law: 



N, d 
T,d 
P 

No,d 

xVdof 



550 ± 40 
361 ±8 
0.41 ± 0.03 

17 ±9 
1161/113 



2±4 
100 ± 28 
0.3 ±0.2 
0.631 ± 2 
1128/296 



40 ±40 
160 ± 20 
0.46 ±0.16 

1±3 
1443/289 



284 ±3 
614 ± 17 
0.00 ±0.01 
0.6 ±0.2 
57/22 



104 ± 1 
8.55 ±5 
0.00 ±0.01 
110 ± 10 
173/25 



12 ±7 
103 ± 11 
0.41 ±0.10 

1±2 
1362/331 




0.001 t 

0.0001 0.0010 0.0100 0.1000 
f. d-' 



100.000 
10.000 
' .000 
0.100 



0.010 
0.001 

0.0001 0.0010 0.0100 0.1000 
f. d-1 




0.0001 0.0010 0.0100 0.1000 1.0000 
f. d-' 



S 0.0100 




O.OCCl I. 

0.0001 C.OOlO 0.0100 0.1000 1.0000 
f. d-i 




0.0001 

0.0001 0.0010 0.0100 0.1000 1.0000 
f. d-' 



10.0000 

^ 1.0000 
in 

D 

^ 0.1000 
0.0100 

fc 

g O.OO'C 



0.0001 \ 

0.0001 0.0010 0.0100 0.1000 1.0000 
f. d-' 



Figure 3: PDS fitting by Lorentzian (upper panels) and modified power law (lower panels) 
models for simulated light curves, models A-C. 
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Figure 4: PDS fitting by Lorentzian (upper panels) and modified power-law (lower panels) 
models for simulated light curves, models D-F. 



4-3. Approximate Solution for Small Variations 

For the linear pseudo- photosphere described by the equation (|5]) , it is easy 
to find the PDS for r in the assumption v = const. For small variations, r, 
r(r) and the resulting luminosity will have identical spectral shapes. 

Let Rq be the radius corresponding to the photospheric radius calculated 
for Mo: 

SktMq 

-Ko - — 

Sirv 

If M[t) is uncorrelated noise with the mean of Mq and dispersion M^Dm <^ 
Mq , the resulting PDS will be identical to the spectrum of the equivalent 5- 
function response: 

2 / — ^1 
rs = -VDrh- 

Here, t > Ro/v, because only visible matter infiuences the opacity. Fourier 
image: 

= fVl>hf (-^ - si(a) + ci(a))) , 
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where integral sine and cosine functions (si and ci) are defined as integrals 
from a = 27rfRo/v to positive infinity. Finally, power spectrum in the linear 
case behaves as: 



The resulting PDS should be peaked, similar to the (|6]) law with p = 
1. The predicted power spectrum is vastly different from the quiet-state 
variability spectrum of LBVs. 

4.4- PDS in the Flaring Case 

The strongly-variable D model has a PDS qualitatively similar to those of 
strongly variable objects like AG Car. In fact, any source exhibiting smooth 
light variations during individual uncorrelated flare events is expected to 
have fiat spectrum at the smallest frequencies and a Brownian-shaped oc 
noise in the opposite limit (/ oo). If one considers a flare of duration tf 
(smooth in its maximum but with rapid rise and decay), its PDS shape in 
the high-frequency limit may be estimated as: 



This approximation works for / ^ In the opposite frequency limit, 
variability is a Poissonian series of short events that implies white noise. The 
spectral slope changes somewhere around / ~ 1/^/, and the exact shape of 
the knee in the PDS depends on the shape of a single flare. The best known 
example is the Lorentzian spectrum for exponential flares. 

Qualitatively, this spectral shape is recognisable in objects like AG Car 
and Romano's star. The steep a ~ 2 slope is probably connected to the 
smooth shape of their flares. Yet the slope remains steep in the lower- 
frequency domain, and the predicted turnover in the power spectrum has 
a frequency a couple of orders of magnitude higher. This may be attributed 
to the existing longer-time variability trends and to correlated behavior of 
outbursts and S Dor cycles. 

Though there are indications for periodicities such as the one-year pe- 
riod in AG Car, a slope holding for two orders of magnitude requires some 
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power- law scaling present. The duration-amplitude relation reported by van 
Genderen et al. ( 1997c ) may be the key. Let us assume that the shape of a 
flare is determined by a single parameter tf (its duration), but its amplitude 
is oc tj. There are indications for p ~ 1. A Poissonian series of such events 
will be characterised by the following power spectrum: 



sin{'Kx)n{x / f)dx 



Here, n{tf) = dN/dtf is distribution of flares per unit duration time. If 
p ~ 1, distribution of the form n{tf) oc tj^ (uniform if one considers Intj 
interval) may explain the observed spectral slopes. 



5. Discussion 

5.1. Different Kinds of Variability? 

One principal question about LBV variability is whether micro- variability 
(with temporal scales less than several months and amplitudes < CP2), S Dor 
cycles, slow luminosity variations and giant eruptions are produced by a 
single physical mechanism. As we show in this article, the basic properties 
of LBV variability at shorter timescales (up to ~ lyr) may be reproduced 
by the toy model described in the previous section. Complex shapes of 
PDS at frequencies / ~ 10^'^ ^ lO^^d^^ (approximately corresponding to the 
duration timescale of an average flare) are possibly connected to transition 
from the wind-smoothed pseudo-photospheric noise at higher frequencies to 
a very similar Brownian noise at lower frequencies, probably connected to 
slow pulsations or some other internal processes. 



Pulsations studied by the hydrodynamical modeling in Fadeyev (2010) are 



similar to the observed micro-variations in periods and amplitudes. They are 
also a good candidate for the process driving mass loss variations. Viewed 
at ~ lyr timescales, irregular pulsations are both fast and uncorrected, that 
justifles the usage of white noise approximation. Longer-timescale variations 
may appear if some instabilities lead to secular evolution of LBV pulsations 



(Buchler, 1993). 



Pseudo-photosphere approach is able to reproduce some of the key points 
of LBV variability such as the steep power spectra. Even if the higher opacity 
of the matter is taken into account, it is short to explain the low- frequency 
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part of the PDSs. This is consistent with the conclusion of de Koter et aL 



(1996) who find that the observed variabihty of LBV stars is too strong to 



be explained solely by velocity and density variations in the wind. 

5.2. Variable Bolometric Luminosity 

Strongest deviations from this simple picture are the giant eruptions and 
"supernova impostors" (see Maund et al. (2006) for review) characterised by 



significant changes in bolometric luminosity (by two magnitudes and more, 



see for example Clark et al. (2009)). Variable bolometric luminosities are im- 
possible to explain in the assumptions of the pseudo-photosphere approach. 

Note however that variable bolometric luminosity may still arise from 
variations in mass loss rate if the characteristic diffusion timescale in the wind 
becomes large enough. Generally, the replenishment timescale ^ is much 
longer than the radiation diffusion timescale that may be roughly estimated 
as tdiff ~ tR/c. This makes the expanding atmosphere practically quasi- 
stationary from the point of view of energy transfer. However, if the total 
optical depth in the wind becomes significantly larger than the inverse wind 
velocity in c units, r > c/f ~ 10'^, the energy output from the pseudo- 
photosphere will be strongly modulated by the envelope evolution. Adiabatic 
losses should be taken into account as well in this case. 

In these energy output modulation, there should be both increases and 
decreases of bolometric luminosity during the outburst. The mean level of 
energy release in the stellar interior is expected to be constant at much larger 
(thermal and nuclear) timescales of > lO^yr. Bolometric luminosity should 
become smaller during giant eruptions because of adiabatic losses close 
to the photon-tiring limit (van Marie et al. , 2009). During strong flares. 



deviations from the mean luminosity level should hold on timescales smaller 
than the survival time of an optically-thick outburst of the mass M expelled 
at the maximal velocity of Vout'- 



1 



out 



Vmit 



KtM, 



out 



4:71 



20 



Vout 



100 kms 



-1 



-1 



M, 



out 



IMo 



1/2 



yr 



(7) 



This estimate naturally arises in a gas that remains hot. Changes in 
opacity and in particular gas recombination produce a variety of similar time 
scales proposed for the photosphere phase or plateau stage durations in super- 



nova light curves (see for example Kasen and Woosley (2009) and references 
therein) . 
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Summarizing, we conclude that LBV eruptions should be accompanied 
by additional energy input inside the hydrostatic radius. 

6. Conclusions 

Variability of LBV stars is generally consistent with a single steep power 
law, but the processes forming its low- and high-frequency parts are probably 
different. The XX-th century light curve of P Cyg does not show any no- 
ticeable flares but is correlated at long timescales up to decades. The overall 
power spectrum is consistent with a p = 1.3 ^ 1.4 power law. For all the 
flaring sources, the slope is higher, p ~ 2. 

Emergence of this steep power-law spectral shape is reasonable in the 
pseudo-photosphere approach. Short-timescale variations in mass loss rate 
are effectively smeared and integrated over the line of sight that leads to 
a Brownian noise with a break at the "replenishment" time of about one 
year. Smaller amplitude variability should produce peaked noise with the 
maximum at several months. The resulting variability lacks two important 
features of the real light curves: it is not correlated at longer timescales (slow 
variability component) and conserves bolometric luminosity. 

Longer-timescale correlations may arise partially from the higher opacity 
of the outflowing matter (that affects the observed luminosity changes in a 
complex way not implemented in our model), variable wind velocity or some 
internal mechanisms. 
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